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We generalize the string method, originally designed for the study of thermally activated rare 
events, to the calculation of quantum tunneling rates. This generalization is based on the analogy 
between quantum mechanics and statistical mechanics in the path-integral formalism. The quan- 
tum string method first locates, in the space of imaginary-time trajectories, the minimal action 
path (MAP) between two minima of the imaginary-time action. From the MAP, the saddle-point 
("bounce") action associated with the exponential barrier penetration probability is obtained and 
the pre-exponential factor (the ratio of determinants) for the tunneling rate evaluated using stochas- 
tic simulation. The quantum string method is implemented to calculate the zero-temperature escape 
rates for the metastable zero- voltage states in the current-biased Josephson tunnel junction model. 
In the regime close to the critical bias current, direct comparison of the numerical and analytical 
results yields good agreement. Our calculations indicate that for the nanojunctions encountered in 
many experiments today, the (absolute) escape rates should be measurable at bias current much 
below the critical current. 
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I. INTRODUCTION 

One of the most important problems in statistical physics involves the rate of decay of a system rendered unstable by 
thermally activated barrier crossing [1-3] and/or quantum barrier tunneling [4-6], and functional integrals represent 
a fundamental tool for studying these transition processes. However, numerical evaluation of the functional integrals 
has always been a challenge. Recently, the string method has been proposed for the numerical evaluation of thermally 
activated rare events [7-11] . This method first locates the most probable transition pathway connecting two metastable 
states in configuration space. The transition rates can then be computed by numerically evaluating the fluctuations 
around the most probable pathway. 

The purpose of this paper is to generalize the string method to the study of quantum metastability caused by 
barrier tunneling. The theory of the decay rate through barrier tunneling has been formulated using the imaginary- 
time functional integral techniques [4] . Essentially, the saddle point of the imaginary-time action is first located and 
the rate of decay is then obtained by evaluating the relevant fluctuations. Within the functional integral formalism, 
the computational task for a quantum field in d-dimensional space is equivalent to that for a classical field in (d + 1)- 
dimensional space [4]. Therefore, the quantum version of the string method can be numerically implemented as its 
original classical version for a higher dimensional system. 
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It should be noted that while in zero-dimensional tunneling problems (a quantum particle is regarded as a zero- 
dimensional quantum field, as in the example presented below) the quantum string method might not offer any special 
advantage than, e.g., the well-known WKB method, in higher dimensional problems or in field theoretical formulations, 
the usual wave- mechanics approach becomes very difficult to implement. It is precisely in such problems that the 
present approach can offer an efficient numerical tool in finding the path of "least action" , and on that basis calculate 
the relevant tunneling rate(s). 

The paper is organized as follows. In Sec. II we outline the string method for the numerical evaluation of thermal 
activation rate and generalize it to the evaluation of quantum tunneling rate. In Sec. Ill we apply the quantum 
string method to study the current-biased Josephson junction [12-15]. This physical model has long been used to 
demonstrate the quantum mechanical behavior of a single macroscopic degree of freedom (the phase difference across 
the tunnel junction) [15]. It has also played an important role in the study of macroscopic quantum tunneling [16]. 
The escape rate of the junction from its zero- voltage state is numerically evaluated at zero temperature in the absence 
of dissipation. For bias currents less than but very close to the critical current, the tilted- washboard potential can 
be approximated by the quadratic-plus-cubic potential, for which the analytic form of the quantum tunneling rate 
has been obtained [17,18]. Our numerical results obtained for this "solvable" model show good agreement with the 
previous analytic results, and thus affirm the validity of the quantum string method. In Sec. IV we conclude the 
paper with a few remarks on the relationship of our results to quantum dissipation. 



II. QUANTUM STRING METHOD 

A. String method for thermally activated rare events 

The string method was originally presented for the numerical study of thermal activation of metastable states [7] . 
Consider a system governed by the energy potential V(q) in the overdampcd regime, where q denotes the generalized 
coordinates {qi}- The minima of V(q) in the configuration space correspond to the metastable and stable states of 
the system. Given and qc as the two minima of V, the most probable fluctuation which can carry the system 
from q^ to qc (or qc to q^) corresponds to the lowest intervening saddle point qe between these two minima, with 
the transition rate given by [1-3] 



T T (A^C)= 1 -^ 

ZTTTj 



|detff(q B ) n " 1/2 
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where r] is the frictional coefficient, H(q) denotes the Hessian of V(q), Xb is the negative eigenvalue of i?(qs), and 
AV = V(qs) — V(qa) is the energy barrier. For the numerical evaluation of r^, we define the minimal energy path 
(MEP) as a smooth curve q*(s) in configuration space. It connects q^ and qc with intrinsic parametrization such as 
arc length s, satisfying 

(VV) ± (cf) = 0, (2) 

where (VV) ± is the component of VV normal to the curve q*(s). Physically, the MEP is the most probable pathway 
for thermally activated transitions between q^ and qc- To numerically locate the MEP in q space, a string q(s, t) (a 
smooth curve with intrinsic parametrization by s) connecting q^ and qc is evolved according to 

(m 1 ^ (3) 

where (dq/dt) 1 - denotes the velocity normal to the string [7]. To enforce the desired parametrization, e.g., equal arc 
length, the string is reparametrized every a few time steps. The stationary solution of Eq. (3) satisfies Eq. (2) which 
defines the MEP. Once the MEP is determined, the intervening saddle point q# is known. The negative eigenvalue 
A b and the corresponding eigenvector can be directly obtained from the MEP, and the ratio of the determinants in 
Eq. (1) can be numerically computed using a stochastic method [9]. 

The above scheme for the calculation of thermal activation rates has many applications, e.g., the condensation 
of a supersaturated vapor, the realignment of a magnetic domain [7,8], and the decay of persistent current in one- 
dimensional superconductor [19,20,11]. All of these transition processes occur when the system undergoes a fluctuation 
that is large enough to initiate the transition. Physically, the thermal activation rate IV becomes practically unmea- 
surable as the temperature is sufficiently low (fcsT <C AV). However, even though thermodynamic fluctuations are 



2 



suppressed at low temperatures, a system can still be rendered unstable by quantum barrier tunneling [4]. The sim- 
plest example is a particle that escapes a potential well: it penetrates a potential barrier and emerges at the escape 
point with zero kinetic energy, after which it propagates classically. In quantum field theory, a classical false vacuum 
is rendered unstable by bubbles of the true vacuum, realized through tunneling. Once these bubbles are sufficiently 
large, they become energetically favorable to grow. 



B. Rate of barrier tunneling 

The theory for the rate of barrier tunneling has been formulated [4] and generalized [5,6] using the imaginary-time 
functional integral techniques. Here we show that the string method can be generalized to be an efficient numerical 
tool (the quantum string method) for the calculation of tunneling rates. In order to make concrete the formulation 
of the approach, below we consider a quantum particle that escapes a potential well through barrier tunneling. This 
zero-dimensional case can be extended to a quantum field, where the classical false vacuum is rendered unstable by 
barrier tunneling [4]. These results show that for a quantum field in d-dimensional space, the computational task may 
be reduced to the calculation of thermodynamic transition rates for a classical field in (d + l)-dimcnsional space. 

Consider a particle of mass m moving in a one-dimensional potential U(q) with two minima qo and q±, one of 
which, qi, is the absolute minimum (see Fig. la). Assume U(qo) = and q\ is to the right of q (qi > q ). The 

normalized harmonic-oscillator ground state i) g (q), centered at qo, is ip g (q) = (mojo/irk) 1 ^ 4 exp [-mwof? — qo) 2 /2h] , 
where uj = yJU"{qo)/m is the frequency locally defined at go- The ground-state energy is hcj /2. These familiar 
results can be derived from the imaginary-time propagator: 

K(q ,q ;T) = (qo\e-" T / h \q ) = J [Dq(r)}e- S / h , (4) 
where T is the imaginary time duration, H is the Hamiltonian of the system, 

S[q(r)} = f' 2 dr{^(^] + U[q(r)} } (5) 

J-T/2 



2 \dr 



is the action, and J[Dq(T)] denotes the integration over functions q(r) satisfying q(—T/2) = q(T/2) = q . Note that 
the imaginary-time (Euclidean) action S^t)] can be obtained from the real-time (Minkowski) action 



A[q(t)] 




dq(t) 12 
dt 



U[q(t)] 



through the formal substitution t —> —it and — z.4[g(i)] — ► S[q(r)]. Thus the equation of motion in imaginary time 
would involve an inverted potential, i.e., U{q) — > —U(q). An expression for K(q ,q ;T) in the limit of T — > oo gives 
both the energy and the wavefunction of the lowest-lying energy eigenstate. In the semiclassical (small U) limit, the 
functional integral for K(q ,q ;T) is dominated by the stationary points of the action, denoted by q(r), that satisfy 
the imaginary-time equation of motion 



SS_ 
~5q~ 



~ TO + u '(«) = °> 



with the boundary condition q(—T/2) = q{T/2) = qo. There are two solutions: one is q(r) = qo, at which the 
Hessian of S, —md 2 + U"{qo), has positive eigenvalues only, and the other is the so-called bounce <7t(r), at which the 
Hessian of S, —md 2 + U"[qb{r)], has a zero eigenvalue plus a negative eigenvalue. The zero eigenvalue comes from 
the time-translation symmetry and the negative eigenvalue makes qb(r) a saddle point of S^t)]. By following the 
bounce %(t) in time, the quantum particle would initially stay at q for a long time, on the order of T, then make a 
brief excursion to the escape point q e (separated from q by a potential barrier, with U(q e ) = U(qo j) in a time of order 
l/wo, and finally return to qo and remain there for another duration of order T (see Fig. lb). This process is called 
a "bounce" . Here q e is the coordinate point at which the quantum tunneling particle leaves the barrier. Physically 
qb{r) characterizes a fluctuation large enough to accomplish the penetration. [For more details, see Appendix A.] 
Using the two stationary points of S in the semiclassical approximation with a proper analytic continuation [4], the 
propagator in Eq. (4) is obtained as 

\Mlo)\ 2 e- E * T/h = (mcooM) 172 e - u ° T ' 2 e- ibaE ' T ' R , (6) 
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which gives l^gilo)] 2 = (mcjo/Trh) 1 ^ 2 , the ground-state energy Tvujq/2, and an imaginary part of the energy ImE g . 
Physically, lmE g (< 0) is responsible for the decay of the metastable "ground state" centered at q 0} with the decay 
rate T Q = -2ImE g /h: 



r<2 = w 



2-nh 



U"(q )\det'H[q b (T) 



detH(qo) 



-1/2 



-s b /h 



(7) 



where H denotes the Hessian of S: H[q(r)] = — md^ + U"[q(T)], Sb = S[qb(r)] is the action associated with the bounce 
<7b(r), and det' indicates that the zero eigenvalue is to be omitted in computing the determinant. 



(a) (b) 




FIG. 1. (a) The potential in one-dimensional space. The unstable ground state is centered at go- After penetrating the barrier, 
the particle emerges at the escape point q e and propagates classically, (b) The bounce solution qb(j) for the imaginary-time 
classical equation of motion. In the inverted potential —(7(g), the particle begins at the top of the hill at go, turns around (i.e., 
bounces) at the classical turning point g e , and returns to the top of the hill. 



C. Numerical implementation of quantum string method 



To numerically evaluate Tq, we generalize the string method to the quantum case. For formal similarity, we write 
the action ^[(/(t)] in Eq. (5) as 5(q), where the vector q represents the coordinates in the g(r)-function space (q 
space). [Computationally, there are always a large but finite number of these coordinates.] We define the minimal 
action path (MAP) as a smooth curve q*(s) in q space. It connects the two minima of S, qo and qi, with intrinsic 
parametrization such as arc length s, satisfying 

(V5) ± (q*)=0, (8) 

where (VS') ± is the component of VS 1 normal to the curve q*(s). Here qo and qi correspond to q{r) — q and 
q(r) = qi, respectively. [A slightly different choice for qi is also possible, corresponding to a q(r) profile with 
q(r) = <?i in most of the r-interval [-T/2, T/2] and q(-T/2) = q(T/2) = q .} The saddle point of S is obtained as the 
point qfc which has the maximum value of S along the MAP. This corresponds to the bounce qb(j~). To numerically 
locate the MAP in q space, a string q(s,t) connecting q and qi is evolved according to Eq. (3) with V replaced by 
S. The stationary solution is the MAP defined by Eq. (8). 

The ratio of determinants in Eq. (7) can be numerically obtained as follows: 

(1) From the MAP q*(s) parametrized by the arc length s, the eigenvector corresponding to the negative 
eigenvalue of the Hessian H(q b ) can be obtained by evaluating dq*(s)/ds at the saddle point q^, followed by a 
normalization. is then computed from A^ = [u^] T W(q(,)u^. 

(2) (2) 

(2) The eigenvector corresponding to the zero eigenvalue X b of the Hessian H(q.b) can be obtained by evaluating 
d T qb(r) from q^,, followed by a normalization. 

(3) The Hessian H(q.b) is modified to give a positive definite matrix 7Y(qh): 

«(<!*) - *(<*) + 2|A«|[u«][u«f + [ui 2 >][uff , (9) 
whose determinant detH(qb) equals | det' H(q.b)\- 

(4) In order to compute the ratio of determinants det7Y(q;,)/det W(qo), a harmonic potential parametrized by a 
(0 < a < 1) is constructed as 
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in the TV-dimensional q space. It can be shown [9] that 



detHK) 
detW(qo) 



exp 



2 [ Q(a)da 
Jo 



where Q(a) is the expectation value 



in which Z(a) is the partition function 

Z(a)= J dqcxp [-ZT(q)]. 
A stochastic process can be generated to measure Q(a) according to Eq. (12) [9]. 



(10) 



(11) 



(12) 



(13) 



III. CURRENT-BIASED JOSEPHSON TUNNEL JUNCTION 
A. Model formulation 

Superconducting devices based on the Josephson effect have been widely used to investigate macroscopic quantum 
tunneling [16]. We consider the resistively and capacitively shunted junction [12,14,15] for which the classical equation 
of motion is 



C 



dt 2 



1 

R 



d<j) 
~dt 



d 



i c h , ih 

COS ( 

2e v 2e 



7i T / N 

= 2e^ 



(14) 



where <j) is the difference in the phases of the order parameters on two sides of the junction, C is the capacitance of 
the junction, R is the resistance of the junction, I c is the Josephson critical current, / is the bias current, and 7jv(t) 
is the fluctuating noise current generated by R. Equation (14) is the same as that of a particle of mass C(h/2e) 2 
moving along the (p axis in an effective potential (the so-called "tilted washboard" potential, see Fig. 2) 



COS0 + — (, 



According to this mechanical analog, for I < I c the zero- voltage state is given by <fio = arcsin(//7 c ), which is a 
minimum of U(<f>). In the classical limit, the escape from the zero- voltage state is induced by the noise current 
which activates the system over the potential barrier [12]. At sufficiently low temperatures, although thermodynamic 
fluctuations are suppressed, the junction can still escape from the zero- voltage state through quantum barrier tunneling 
[14,15]. 



CD -2 

^ -4 - 




10 



20 



4> 



FIG. 2. The tilted washboard potential U{4>) for I/I c = 0.3. Here the unstable ground state is centered at </>o, <f> e is the 
escape point, and (f>\ — <f>o + 2tt is the next minimum. 
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We consider the zero-temperature behavior of the junction in the low-damping limit ((2e/ c /?iC) 1 / 2 [l — 
(I// c ) 2 ] 1 / 4 i?C 3> 1). The state of the junction is represented by the wave function ^(4>, t), governed by the Schrodinger 
equation ifid^i jdt — in which the Hamiltonian is of the form 



We introduce three parameters: the dimensionless bias current x = I/I c , the plasma frequency uj v = y/2eI c /hC, and 
I c /2eu) p = \J Ej I Ec where Ej = I c fi/2e is the Josephson coupling energy and Eq = (2e) 2 /C is the charging energy. 
The imaginary-time action 



S[4>(r)} = J 



dT 



:C 



lb- 



+ U{4>) 



(15) 



can be written as S[0(t)] = Ej / 'EcS[</>(t)], where r = uj p t is a dimensionless time variable and S is the dimen- 
sionless action 



S[0(f)] 



df 



1 

2 V df 



+ (— cos< 



(16) 



B. Numerical results 

Based on the actions (15) and (16), the tunneling rate is given by 



T Q = UJp^COS O (j^ 




- — 1/2 

cos0 o |det H[4>b(f)]\ 
del H(ou) 




exp -J-ET-St , (17) 



according to the general expression (7). Here uj p \/ cos O is the frequency locally defined at <j>o, Sb = <S[0{,(f)] is 
the dimensionless action of the bounce 0b(f) which satisfies 0b(±oo) = 4>o, and 7i is the Hessian of 5[^(f)], given 
by H[4>(f)] = —d 2 /d? 2 + cos[^(f)]. Numerical calculations based on the action in Eq. (16) have been carried out 
to evaluate the bounce <j>b(f), the bounce action Sb, and the ratio of determinants in the general expression (17). 
Note that these dimensionless properties are uniquely determined by the parameter x. Once they are evaluated, 
the tunneling rate Tq can be readily obtained using the other two parameters oj p and y/ Ej / Ec ■ These parameters 
should be easily measurable experimentally, since they are directly determined from the capacitance of the junction, 
the critical current, and the bias current. 

Numerical calculations have been carried out according to the following procedure. 

(i) We first locate the MAP in the </>(f)-function space. In the calculation, 0(f) is represented by a column vector <fi 
of N = 200 entries, with the f-interval [-T/2, T/2] discretized by a uniform mesh of N points. We use T = 20, large 
enough for the computation of zero-temperature properties. [Here T = ftwp/fcsT, where ks is the Boltzmann constant 
and T the temperature, and T ^> 1 means hui p ^> ksT.] The string <fi(s) connecting 0(0) = <p and 0(1) = (fr-^ is 
discretized by M = 200 points in the 4> space. As to the two fixed ends of the string, cf> corresponds to 0(f) = O 
and 4>\ to a 0(f) profile with 0(f) = 0i in most of the f-interval and 0(— T /2) = 0(T/2) = O . Here O = arcsinx 
and 0i = 0o + 2ir are two neighboring minima of — (cos0 + x0) (see Fig. 2). Note that <p 1 is obtained as a local 
minimum of S in Eq. (16). The string evolution is generated by 

(d0/dt)- L = -(VS) ± (0) ) 

with the initial string taken from a linear interpolation between 4> and <f> 1 . The MAP <fi*(s) is reached by the evolving 
string <p(s,t) as its stationary solution defined by (V5) (</>*) = 0, with 0*(O) = 4> and 0*(1) = 4>i- The bounce 
0b (f) is obtained from the vector (f> b which yields the maximum value of S along the MAP. In Fig. 3 a sequence of the 
0(f) profiles along the MAP is shown for x = 0.1, and in Fig. 4 the bounce profile 0&(f) is shown for a few selected 
values of x. In Fig. 5 the variation of the action S along the MAP is shown for a few selected values of x, and in Fig. 
6 the bounce action Sb is plotted as a function of x. 



6 





FIG. 4. The bounce <f>b (r), plotted for a few selected values of x. 
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0.2 0.4 0.6 0.8 



S 

FIG. 5. The dimensionless action along a segment of the MAP starting from 0(r) = (j>o and ending at 4>b(j). Here 5* is 
plotted as a function of the arc length s in the ^(f)-function space for a few selected values of x. The profile of <f>(f) = 0o is 
taken as the reference point s — at which S has been set to be zero by a constant shift of the potential. 
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0.2 0.4 0.6 0.8 1 
X 

FIG. 6. The bounce action Sb plotted as a function of x. 

(ii) We then modify the Hessian H[^>b(f)], represented by the N x N matrix H(4> b ), to give a positive definite 
matrix H((fi b ), given by 

H(0 6 ) - H& b ) + 2|A[ 1) |[u[ 1) ][u[ 1) ] T + cos0 o [ui 2) ][u[ 2) ] T . 

Here is the eigenvector corresponding to the negative eigenvalue A^ of the Hessian 7i(4> h ), is the eigenvector 
corresponding to the zero eigenvalue X b of the same Hessian, and det7i(</> 6 ) = cos 4>q\ det' 7i(4> b )\. Note that 
A^, and can be readily obtained once the MAP is determined, as outlined in Sec. II C. In Fig. 7 the unstable 
direction along the MAP at the saddle point is shown for a few selected values of x. 




FIG. 7. The eigenfunction u£ (f) of the Hessian Tl[<f>b(f)] with the negative eigenvalue, plotted for a few selected values of a;. 
The curves are displaced vertically for clarity, whereas the original ones all start and end at 0. Each eigenfunction represents the 
unstable direction at the saddle point along a particular MAP in the </>(f)-function space. It is noted that uj; 1 ' (r) obtained for 
x = 0.1 is qualitatively different from that for x — 0.9. For x close to 0, the bounce <j>b{j) has the height <f> e (the escape point) 
close to the next (lower) minimum 0i, and the growth of the 4>(f) bubble along the MAP is characterized by the movement of 
the two (left and right) domain walls. Consequently, the unstable direction u^(f) shows two peaks. On the other hand, for 
x > 0.3, the bounce 4>b{f) has the height <j> e far from the next minimum cj>\, and the growth of the <j>{j) bubble along the MAP 
is characterized by overall dilation. Consequently, the unstable direction u b (f ) displays one peak only. 



(iii) We calculate the ratio of determinants 



_ cos 0o I det' H [(fa (t)] I _ cos <f> \ det' H{4> b )\ 
7 ~ detH((/) Q ) ~ detH(0 o ) ' 1 ' 
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where 7i{<p ) is the N x N matrix representation of TC((f)o). This is done by evaluating det T-i(<p b ) / det 7i.((f) ) according 
to the stochastic method outlined in Sec. II C. Numerical results of this part are shown in Figs. 8 and 9, for the 
stochastically measured Q{a) in Eqs. (11) and (12) and the ratio of determinants 7 in Eq. (18). From Sb and 7, the 




0.2 0.4 „ 0.6 0.8 1 



FIG. 8. Stochastically measured Q(a), plotted for a few selected values of x. 
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FIG. 9. The determinant ratio 7 in Eq. (18) plotted as a function of x. 

Using the numerical results for Sb and 7, the mean escape rate out of the zero- voltage state can be readily obtained 
from Eq. (17), once the values of w p and Ej / Ec are given. From I c = 9.489 fiA and C = 6.35 pF reported in an 
early experiment [15], we have tu p — 67.4 GHz and \J Ej /Ec = 440. The largeness of WEj/Ec implies that quantum 
tunneling becomes observable only if x — ► 1 and Sb — > 0. For the experiment reported in Ref. [15], x ~ 0.99 at 
which Sb ~ 0.037, e~v ' E -> '/ IE c.s b ^ 10~ 7 , and the escape rate is approximately 2.7 x 10 4 sec -1 . In a recent experiment 
on quantum superposition of macroscopic persistent-current states, a superconducting loop is constructed with three 
Josephson junctions [21]. We find that the junction parameters in that experiment allow quantum tunneling to be 
observable in a range of x much wider than that in Ref. [15]. From I c = 570 nA and C — 2.6 fF [21], we have 
ujp = 816 GHz and y/ Ej / Ec — 2.18. The smallness of y Ej /Ec then allows Sb and hence x to vary in a wide 

range. For x decreasing from 0.8 to 0.2, Sb roughly increases from 2 to 10, and consequently, e~^ Ej ^ Ec ^ b changes 
from ~ 10 -2 to ~ 10~ 10 . Using Eq. (17) with the numerical results for Sb and 7, we obtain the escape rates 
Tq = 1.2 x 10 11 sec -1 , 7.4 x 10 7 sec" 1 , and 1.5 x 10 3 sec -1 , for x = 0.8, 0.5, and 0.2, respectively. Note that typically, 
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the measured escape rates are in the range from ~ 10 1 to ~ 10 6 sec" 1 . Therefore, our numerical results indicate 
that the absolute escape rates for today's nanojunctions should be measurable at bias current much below the critical 
current. This is because the junction capacitance has been significantly reduced and thus the action scale y/Ej/ Eq 

can be made small enough to allow a relatively large dimensionless action St, in the exponential factor e~V- Ej /- Ec- ^ 
[22] . In this regard a numerical scheme as presented in this paper is essential to the evaluation of the absolute escape 
rates. 




FIG. 10. The dimensionless prefactor y'cos 4>o\ — 
the large variation of a factor > 10 as a function of x. 



^o\det' H[<t> b {f)]\ 
det"W(0 o ) 



-1/2 



in Eq. (17), plotted as a function of x. Note 



C. Quadratic-plus-cubic potential 

For bias currents less than but very close to I c (x close to 1), the potential barrier AU to be penetrated is low, the 
distance between the minimum 0o and the escape point e is small, and hence the potential U(cj>) in the classically 
forbidden region can be approximated by the quadratic-plus-cubic potential. That is, the dimensionless potential 
— cos <fi — xcf) in the action (16) can be expressed in a Taylor expansion form 

u{4>) = ~ cos - x(f) fa (- cos O ^ £0o) + i cos O (0 - <j) ) 2 ( 1 - — — -p J , (19) 

2 \ <pe — <pa J 

where e — 0o = 3cot0o approaches zero as x — > 1 and (f>o — * 7r/2. From Eq. (19), the dimensionless potential 
barrier can be easily found to be Au = 2cos</>o(^ e — </>o) 2 /27 and the dimensionless bounce action to be Sb = 
8Vcos (f>o{(f>e — 0o) 2 /15- In addition, the Hessian of S{4>(?)] becomes H [0(f)] = — 9 2 /9f 2 +cos0 o [1^3(0— 0o)/(0e — 0o)]> 
for which the analytic result 

cos 0o I det'W[0b(f )]| _ 1 
detH(0 o ) _ 60 

has been obtained for the determinant ratio 7 in Eq. (18) [17,18]. These exact results for the quadratic-plus-cubic 
potential allow an analytic form for the tunneling rate Tq in Eq. (7). Numerical results for x = 0.9 in Figs. 6 and 

A 

VCOS 00 (0e ~ 00 ) 2 

though e — 0o = 1.45 is still large. 

In order to demonstrate the validity and precision of the quantum string method, numerical calculations have been 
carried out to reproduce the bounce action and determinant ratio for the quadratic-plus-cubic potential, an important 
model potential for the study of quantum metastability [18]. For simplicity, we work with the scaled action functional 



9 show that — = 0.469, approaching 8/15, and the determinant ratio 7 = 0.0162, approaching 1/60, 



S[q(r)} = / dr 



(20) 
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The potential g 2 (1 — g) /2 in action (20) has qo — as the metastable minimum and q e = 1 as the escape point. For 
computational purpose, this potential is distorted in the region of q 3> q e to generate another (lower) minimum at 
Qi (^ Qe)- Numerically, the two potential minima qo and q\ are used to fix the ends of the evolving string in the 
g(r)-function space. The stationary solution for the string evolution equation is the MAP from which the saddle point 
of S[q(r)], i.e., the bounce Qb(r), can be obtained. Since the potential profile is untouched in the classically forbidden 
region (go < q < qe), the bounce so obtained is not affected by the potential distortion far away. The bounce action 
Sb = S[qb(r)] is obtained to be 0.5337, very close to the exact result 8/15. The determinant ratio 

\det' H[q b (r)}\ = |dct'[-d 2 /dT 2 + l-3 gb (r)] | / 1 \ 
detH(g ) dct(-d 2 /<9r 2 + l) \ J 

is obtained to be 0.0142, close to 1/60. These results are obtained for N = 200, M = 200, and the total imaginary 
time duration T = 20. Better agreement with the exact results can certainly be achieved by using longer imaginary 
time duration, vector space of higher dimensionality, and finer resolution in discretizing the string. 



IV. CONCLUDING REMARKS 



Quantum tunneling in macroscopic systems is intimately related to the important issue of quantum dissipation, 
which arises from the coupling to environmental variables. This coupling can modify the tunneling itself, as many 
prior works have shown [5,17,23,24]. However, it should be noted that regardless of the particulars in the quantum 
dissipation model, the net result is to decrease the escape rate. Hence the rate calculated for nondissipative quantum 
tunneling may be regarded as an upper bound to the rate(s) with nonzero dissipation. 

In this regard it should also be noted that as a tool for the numerical evaluation of tunneling rate in the path 
integral formalism, the quantum string method is directly generalizable to field theoretic problems, requiring only 
additional computational resources. 
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APPENDIX A: THE BOUNCE 



The bounce <7&(t) is a solution of the imaginary-time classical equation of motion: 

m— = U (q), 

subject to the boundary condition q(—T/2) = q(T/2) = qo for T — > oo. The qualitative behavior of qb{r) is suggested 
by the analogy with the equation of motion for a particle of mass m in the inverted potential — [/(g), in which q now 
corresponds to the top of the hill and q e to the classical turning point. The particle would spend most of its time at 
go (due to zero velocity), but, in the course of an arbitrarily long interval of time, it would make a brief excursion to 
the point q e and then return to go (see Fig. lb). Note that 

^( d 4) 2 -U[Q(r)] = -U {qo ) = 



2 \dr_ 

is a constant of motion for gb(r). This means dqb/dr vanishes at go and q e . 

The bounce gb(r) shown in Fig. lb is centered at r c = along the r axis. Because of the time-translation invariance 
of the action, the bounce solutions are also given by qb(r — r c ), where r c is an arbitrary center of the bounce. This 
symmetry property leads to a zero eigenvalue for the Hessian of S at gb(r), — ra<9 2 + {/"[<#, (t)]. The corresponding 
eigenfunction is given by 

u (2 Ur)- 

Ub [T) -ys b dT> 
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where ^Jm/Sb is the normalization factor derived from the action of the bounce, 




and the constant of motion 

!(t) 2 -"M^°- 

Note that dqb/dr has a zero at the center of the bounce. Therefore, u b (r) has a node and can not be the eigenfunction 

with the lowest eigenvalue: there must be a nodeless eigenfunction, u^(t), with a negative eigenvalue. This implies 
that the bounce is not a minimum of the action but a saddle point. The negative eigenvalue requires a proper 
analytical continuation in evaluating the functional integral in (4). This leads to a complex energy in (6). 
Using the constant of motion, the bounce action can be reduced to the form 

S b = J ' 1 m (^j dr = lJ* C y/2mU{q)dq. 

It is seen that e ^ Sb ^ 2h is the familiar WKB exponential factor for the amplitude of tunneling wave. Accordingly, 
e -s b /h - g exponential factor for the current density of the tunneling wave, which is directly related to the rate of 
decay. This is reflected in Eq. (7). We want to remark that for one-dimensional quantum mechanics, the tunneling 
rate (7) derived from functional integral agrees with that obtained by standard WKB method of wave mechanics [4]. 
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